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Abstract 

The sokition to the nonhnear output regulation problem requires one to solve a 
first order PDE, known as the Francis-Byrnes-Isidori (FBI) equations. In this paper 
we propose a method to compute approximate solutions to the FBI equations when the 
zero dynamics of the plant are hyperbolic and the exosystem is two-dimensional. With 
our method we are able to produce approximations that converge uniformly to the true 
solution. Our method relies on the periodic nature of two-dimensional analytic center 
manifolds. 

1 Introduction 

Consider the control system 

X = f{x, u, w) 

lit = s{w) (1) 
y = h{x, u, w) 

where x G is the state variable, u G is the control variable, G M"^ is an exogenous 
variable, and y eW> is the output variable. The maps / : M" x x ^ M", s : ^ M"? 
and h -.W^x x M"? — )■ Rp are all assumed to be sufficiently smooth and satisfy /(O, 0, 0) = 0, 
s(0) = 0, and /i(0,0,0) = 0. The variable w represents a disturbance and/or a reference 
signal, and its dynamics are commonly referred to as the exosystem. The state feedback 
regulator problem [8] is to find a state feedback control u = a{x,w), with a(0,0) = 0, such 
that the equilibrium a; = of the dynamical system 

X = f{x, a{x, 0), 0) 

is exponentially stable, and such that for each sufficiently small initial condition (xq, Wq) the 
solution of dl]) with u = a{x,w) satisfies 

lim y{t) = 0. 

A characterization of the state feedback regulator problem for linear systems was given 
by Francis [3] and later generalized to nonlinear systems by Isidori and Byrnes [8]. As 
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shown in [8], the solvabihty of the regulator problem can be reduced to the solvability of 
a system of partial differential equations (PDEs), which in the linear case reduce to the 
Sylvester type equation obtained by Francis. For this reason, we refer to these equations as 
the Francis-Byrnes-Isidori (FBI) PDEs. For completeness, we state the main result of [8] 
(for the definition of Poission stability used below see Remark II. ip . 

Theorem 1.1 (Isidori-Byrnes) . Assume that the equilibrium w = of the exosystem is 
Lyapunov stable and there is a neighborhood of w = in which every point is Poisson stable. 
Assume further that the pair 



is stabilizable. Then the state feedback regulator problem is solvable if and only if there exists 
(k > 2) mappings n : n ^ W, with 7t{0) = 0, and k : n ^ W^, with k(0) = 0, both 
defined in a neighborhood f2 C M'^ of w = 0, and satisfying 



Given a solution pair (tt, k) to the FBI equations ([3]), a state feedback solving the regulator 
problem is given by 



where K G x R" is any feedback matrix rendering the pair ([2]) asymptotically stable. In 
general, solutions to the FBI equations, being singular quasilinear PDEs with constraints, 
may not exist. However, for a class of control- affine systems, it is shown in [8] that the 
solvability of the FBI equations is a property of the zero dynamics of ([T]). Roughly speaking, 
if the zero dynamics of ([T]) has a hyperbolic equilibrium at the origin then a solution to 
the FBI equations exists by the center manifold theorem [2]. It is known, however, that 
center manifolds suffer from a number of subtle properties associated with uniqueness and 
differentiability [I2]. Despite these difficulties, a C°° dynamical system possess a C*^ center 
manifold for each k > 1, and moreover, it is possible to obtain approximate solutions of 
arbitrarily high-order via Taylor series [2j. In this respect, Huang and Rugh [B] and Krener 
^ provide a method to compute approximate solutions to the FBI equations via Taylor 
polynomials, yielding approximate output regulation. A shortcoming of this approach is that 
the domain on which the series approximation yields satisfactory results is not guaranteed 
to enlarge significantly by computing higher order approximations. This can be a serious 
drawback as the number of monomials in q variables of degree d is (''^j"^) , a number growing 
rapidly in d. Moreover, polynomial approximations to (tt, k) can lead to destabilizing effects 
when the order of the approximation increases. 




(2) 



dw 



{w)s{w) = f {71(10), k{w),w) 
= h{Ti{w), k{w), w). 



(3) 



a{x, w) = k{w) + K(x — -K^w)) 
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In this paper, we present a method to compute solutions to the FBI equations for the 
class of real analytic SISO control-affine systems 

X = f{x) + g{x)u 

w = s{w) (4) 
y = h{x) + p{w) 

where / : M" ^ M", ^ : ^ M", s : M« ^ W, h : W ^ M, and p : M« ^ M are 
real analytic about the origin. Furthermore, we will restrict our considerations to two- 
dimensional exosystems, i.e., q = 2, whose linear part contains non-zero eigenvalues. Our 
method is based on the existence and uniqueness results for two-dimensional analytic center 
manifolds in yj and on the high-order patchy method of Navasca and Krener [lOj. A key 
strength of our approximation method, which relies on the periodic nature of the solution to 
the FBI equations for (jl]), is the reduction of the computational effort inherent in a direct 
Taylor polynomial approximation. 

The organization of this paper is as follows. In Section [2] we briefly summarize the 
key insight provided in [S] on how the solvability of the FBI equations can be reduced to 
the problem of solving a center manifold equation provided the zero dynamics of (jlj) are 
hyperbolic. With this simplification, we show how the standard stability assumptions on 
the exosystem lead to a direct application of the results in [Ij to deduce uniqueness of 
solutions to the FBI equations of (j4]). In Section [3l we describe our method to compute a 
high-order piecewise smooth approximation to the solution of the FBI equations and prove 
that a sequence of approximations generated by our method converges uniformly to the true 
solution. Finally, in Section H] we illustrate our method on examples and then make some 
concluding remarks. 

Remark 1.1. Henceforth, it will be implicitly assumed that the exosystem has an equilib- 
rium w = that is Lyapunov stable and that there is a neighborhood of = in which 
every point is Poisson stable. We will refer to this type of stability as neutral stability. By 
Poisson stability we mean the following. An initial condition xq of the dynamical system 
X = f{x) is Poisson stable if the flow ^{{xq) of the vector field / is defined for alH G M and 
for each neighbourhood U of xq and for each real number T > 0, there exists a time ti > T 
such that ^{^ (xo) G U and a time t2 < —T such that ^{^ (xq) G U. 

2 Real analytic and periodic solutions to the FBI equa- 
tions 

As shown in [8], a key simplification in the problem of solving the FBI equations for a system 
of the form (jl]) consists in reducing it to the problem of solving a center manifold equation 
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for the zero dynamics of (jl]). Following [8] and using the now standard notation in 
assume that the triplet {/, g, h} has relative degree l<r<'n,atx = and let {z, ^) denote 
the standard normal coordinates, where ^ = {h{x), Lfh{x), . . . , Uj~^h{x)) and z is such that 
LgZ = 0. In the {z,C,) coordinates, (jlj) takes the form 

^1 = ^2, ■ ■ ■ , Cr-l = 

ir = b{z,0 + aiz,Ou (5) 
ib = s{w) 

The zero dynamics of (j4]) are given by the dynamical system 

i = /o(^,0). (6) 

Define functions (y^j : — )■ M by i^iiw) = —U^^p{w), 1 < i < r, set ^{w) = {ipi{w), . . . , ipr{w)), 
and let 

Ue{x, W) ■ 



L)h{x) + Llp{w] 



LgUph{x) 

Then it is straightforward to verify that if (j) satisfies the PDE 

— {w)s{w) = /o(0(w), (/?(«;)) (7) 

then tt{w) := {(l){w),(p{w)) and k{w) := Ue{7r{w),w) constitute a solution pair to the FBI 
equations of ([5]). If the origin of is hyperbolic, then ([7j) is the equation that is satisfied 
by any center manifold {{z,w) : z = (piw)} of the dynamical system 

z = fo{z,^{w)) 

[o) 

w = s{w). 

Hence, in the hyperbolic case, the problem of solving the FBI equations associated to the 
original system (jl]) is reduced to solving the center manifold equation associated to 
Although this simplification is significant, solutions to center manifolds suffer from subtleties 
associated with uniqueness and differentiability [T5]. For example, it is known that an 
analytic dynamical system does not generally posses an analytic center manifold, thereby 
forcing one to seek a center manifold solution that is only {k = 2,3,...) and thus not 
necessarily unique. As an example, the polynomial dynamical system 

z = —z + wl + wl 

Wi = -W2 - lwi{wl + wl) 
W2 = Wi- ^W2{wl + wj) 
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which is of the form ([8]), has the property that each center manifold has Taylor series 

oo 

1=1 

which has vanishing radius of convergence. Despite these difficulties, a special case for which 
sharp uniqueness and differentiability results exist is for two-dimensional center manifolds, 
and is given by the following theorem due to Aulbach [1]. 

Theorem 2.1 (Aulbach). Consider the ordinary differential equation 

z = Bz + Z{wi, W2, z) 

Wi = -W2 + P{Wi,W2, z) (9) 
W2=Wi + Q{Wi, W2, z) 

where Wi,W2 G M, -2 G M", and P,Q, and Z are real analytic functions about the origin 
and have Taylor series beginning with quadratic terms. Suppose that the matrix B has no 
eigenvalues on the imaginary axis. If the local center manifold dynamics of ([9]) are Lyapunov 
stable and non- attractive then ^ has a uniquely determined local center manifold which is 
analytic and generated by a family of periodic solutions. 

Aulbach's result has a direct application to the output regulation problem, as given by 
the following theorem. 

Theorem 2.2. Suppose that in (jl]) the exosystem is two-dimensional and ^(0) has non- 
zero eigenvalues. Suppose that f,g,h and p are real analytic mappings about x = and 
w = 0, respectively, and that the triple {f,g, h} has a well-defined relative degree 1 < r < n 
at X = 0. If the zero dynamics of (jl]) are hyperbolic, then there exist unique and real analytic 
mappings (tt, k) solving the associated FBI equations of (jl]). 

Proof. By assumption and neutral stability of the exosystem, the eigenvalues of the exosys- 
tem are non-zero and purely imaginary. Indeed, if the eigenvalues were not purely imaginary 
then w = would necessarily be either a repelling or an attractive equilibrium, contradicting 
the assumption of neutral stability. Now since B := ^(0,0) contains eigenvalues off the 
imaginary axis, there exists an analytic coordiante change [H] about the origin such that 
(jH]) takes the form 

Z = Bz + Z{wi, W2, z) 
Wi = -W2 + P{Wi,W2) (10) 
W2=Wi-\- Q{Wi, W2), 

where P, Q and Z are analytic at the origin and have Taylor series beginning with quadratic 
terms. From (IIUI) . we can observe that the dynamics of any center manifold of (llUp are 
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equivalent to the exosystem dynamics, which by assumption are Lyapunov stable and non- 
attractive. Aulbach's theorem completes the proof. ■ 

Remark 2.1. Theorem 12.21 actually holds for more general MIMO control-affine systems 
with m = p. In [5j , it is shown that if the composite control-affine system 

m 

X = f{x, w) + '^ gi{x, w)ui 

i=l 

w = s{w) 
y = h{x, w) 

has a well-defined relative degree at {x,w) = (0,0), then the associated FBI equations are 
solvable if the zero dynamics of the composite system are hyperbolic. In this case, the FBI 
equations reduce to a center manifold equation of the form ([7]) so that Aulbach's theorem 
can be applied when the exosystem is two-dimensional. 

Example 2.1. The dynamics of a cart and inverted pendulum system can be written in the 
form 

Xi = X2 
X2 = U 

(11) 

Xs = X4 ^ ' 

±4 = 1 sin(x3) - J cos(x3)w 

where Xi is the position of the cart, x^ is the angle the pendulum makes with the vertical, g 
is the acceleration due to gravity, £ is the length of the rod, and u is the control force. With 
h{x) = Xi, the system has relative degree r = 2 at x = 0, and therefore (^1,^2) = C(^) = 
{h{x), Lfh{x)) = (xi,X2). With (21,-22) = z{x) = (X3, X4 + cos(x3)), the zero dynamics are 
given by 

Zl = Z2 

Z2 = -sm(zi), 

whose linearization has eigenvalues ±A/f- Hence, with system output y = Xi + p{w) {p real 
analytic) and a two-dimensional real analytic exosystem whose linearization has non-zero 
eigenvalues, there exists a unique and real analytic solution to the associated FBI equations 
of the cart and inverted pendulum system ( fTTj) . □ 

3 Computation of the center manifold 

In this section we outline a method to compute the solution to the FBI equations in the case 
of two-dimensional exosystem and real analytic data. As described in the previous section, 
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for the nonhnear control systems in consideration, the solvability of the FBI equations can 
be reduced to solving a center manifold equation for a dynamical system of the form 

Z = Bz + Z{wi, W2, z) 
Wi = -W2 + P{Wi,W2) (12) 
W2=Wi + Q{Wi, W2), 

where w = (101,^2) G M^, z G M", Z, P, and Q are real analytic mappings, and the 
eigenvalues of B have non-zero real parts. We will therefore limit our considerations to 
solving the center manifold equation for f[T^ . It will be assumed that the ty-dynamics have 
w = as a Lyapunov stable and non-attractive equilibrium. By Theorem 12. H there exists a 
unique analytic mapping 4>{wi, W2), defined locally about w = 0, solving the center manifold 
PDE associated to ffT2|) . 

Our method is best described on the representation of (fT2l) in polar coordinates. Hence, 
we apply the tranformation {wi,W2,z) = {r cos6,r sinO, z) to f[T^ yielding a system of the 
form 

r = rR{6, r) 

^= 1 + 6(6/, r) (13) 
i = Bz + Z{e,r, z) 

where R,Q,Z are analytic functions converging for each 9 G [0,27r] and |r| < a, ||2;|| < a, 
where a > is a positive constant. Define f{9, r, z) = Bz + Z{9, r, z). The center manifold 
PDE for ([13]) is 

Mr,^(^,r)) = ^[l + e(^,r)] + ^rR{9,r) (14) 

for the unknown analytic mapping '0(6', r) (= 0(r cos 6, r sin 6)). The mapping ip has a power 
series representation 

00 

'i|J{9,r) = Y,^^{0y 
1=1 

converging in a cylinder of the form 9 G [0,27r], |r| < e, and with 27r-periodic coefficients 
Gi{9) [I]. By eliminating the time variable t, f[T^ can be reduced to 

fir 

- = rR{9,r) (15a) 
= Bz + Zi9,r,z). (15b) 

a9 

Define f{9, r, z) = Bz + Z{9, r, z). From it follows that 
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We now give a brief sketch of our method. Let r{9) be a solution to fllSap and define the 
mapping 

for ^ G [0, 27t] and \a\ small. We note that, with ; r = r{9) substituted into the RHS of (fTSbj) . 
the curve \l/(6', 0) = ■ip{9, r{9)) is the solution to fllSbp with initial condition z{Q) = ip{0, ^{^))- 
For \a\ sufficiently small, we have a power series representation 

m^) = no,o) + Y,^{9,o)^ (17) 

i=l 

converging for all 9 G [0,27r] and having 27r-periodic coefficients |^(6',0). In fact, it is easy 
to see that 

a^(^'0)=a^(^'-W)- (18) 

By construction, the mapping \1/ is a perturbation of 'i/j{9,r{9)) in the radial direction, the 
amount of perturbation given by the parameter a. Our method is based on computing the 
Taylor series approximation 



1=1 



and using it to build the center manifold along r{9) in the radial direction. Having followed 
along a small annular region, say of the form 

{i9,r):0<9< 2n, r{9) <r < r{9) + e}, 

we compute a new radial curve 9 i— )■ f{9) with initial condition f(0) = r(0) + e, compute the 
new corresponding Taylor series approximation \E'^, and then continue building the center 
manifold by following along the annular region 

{i9,r):0<9< 2n, f{9) <r < f{9) + e}. 

This process is repeated and the annular regions, along with the corresponding approxima- 
tions, are patched together to form a piecewise smooth approximation to the true solution 
^. 

Remark 3.1. To compute the Taylor series approximations it is necessary to compute 
the ^-dependent coefficients appearing in (fT7|) . which can be done in the following way. From 
the definition of a direct computations gives 
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which when combined with ( !T6|) yields 

^ = f{e,r{e) + a,^{d,a)). (19) 

Using (fT9l) . we can now write a down hnear inhomogeneous ODE for the coefficient |^(^, 0). 
Indeed, differentiating ( fT9l) with respect to a, and interchanging the order of differentiation, 
yields 

d f \ df d^i df 



and therefore 



5 ( \ df 



where the matrix A{9) = ^{0, r{9), "^{9, 0)). In general, it can be verified by induction that 

m (a^rC. 0)) = m^io, 0) + f, [e, me, 0). ^(m). .... g-r^-. 0)) (20) 

for some mappings Fj, i > 2. 

With the previous constructions in mind, we are now ready to describe an algorithm for 
computing the solution -0 to the center manifold equation ( |T6l) . 

1. Let > 1 be a fixed positive integer and let 



N 



i=l 

that is, -01^ is simply the A^th order Taylor approximation of -0 in r. To compute i/jq , one 
can use the method in to generate a A^th order Taylor polynomial approximation of 
(f){wi, W2), say (f)'^{wi, W2), and then simply i/jq {9, r) = 0^(r cos^, r sin^). Set \I'o = "0 
and set r_i{9) = for 9 E M.. The initial approximation ■0^ will be accepted in an 
annular region of the form 

{{e,r) :0<^<27r,0<r< ro{9)} 

where ro(6') is the solution to f llSap with some prescribed initial condition ro(0) = eo > 
0. To compute accurate numerical solutions to tq, we solve a BVP using fllSap with 
boundary conditions r(0) = r(27r) = eo and constant initial guess eo on [0,27r]. 

2. Define cr) = tp{9, rQ{9) + a). From f[T71) . \E'i can be approximated by the truncated 
series 

^n0,-) = M0,0) + f2^{9,0)^ 

i=l 
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for \a\ small. To obtain accurate numerical solutions to the coefficients ^^(6',0), we 
solve EVPs using the ODEs ^ with boundary conditions ^(0,0) = ^(27r,0) 
and initial guesses 

Similarly, to compute 0) = ip{0, ro{9)) we solve a EVP using fllSbl) with boundary 
conditions z{0) = z{2tt) and initial guess 

M/i(^,O)^<(0,ro(^)). 

Having computed ^i{6, 0), ^^{d, 0), . . . , ^^(6*, 0), we obtain an approximation ipii^O, r] 
to il'{0,r) defined by 

which is accepted in the region 

{{e,r):0<e< 27r, ro{9) <r< ro{9) + (21) 

for some desired ei > 0. In this way, we have extended our original approximation i(jq 
of ip to the domain (12T!) . Our running approximation of i/j is given by 



>o(^,r), 0<r<ro{9), 

ro{e) < r < ro{e) + e. 



for ^€[0,2 



vr 



3. We now proceed to augment to our running approximation a mapping ijj2, that will be 
defined on an an annular region surrounding the domain of in the following way. 
We first compute the solution ri{9) to fllSap with initial condition ri(0) = ro(0) + ei. 
As in Step 2, this is done by solving a EVP using (I15al) with boundary conditions 
r(0) = r(27r) = ro(0) + ei and taking the curve ro(-) + ei as an initial guess to ri. 
Here we note that, to avoid overlapping domains of definition between i/ji and 'ijj2, the 
domain fl^ of ipi is redefined to be 

{{9,r) ■.0<9< 27r,ro(9) <r < ri{9)}. 

4. We now repeat Step 3 with ri{9) and build an approximation to "^2(9, o") = 1^(9, ri{9) + 
cr) of the form 

^no,a) = Mo,o) + f2^{9,o)^, 

i=l 
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for < cr < £2 and £2 > sufficiently small. The coefficients ^^{9, 0) are computed by 
solving EVPs using the ODEs ([2DD with boundary conditions f^(0,0) = ^(27r,0) 
and initial guesses 

That is, we use the previously computed coefficients as initial guesses for the current 
coefficients. Similarly, to compute "^2(9,0) = 'i/j{9,ri{9)) we solve a EVP using ( llSbl) 
with boundary conditions z{0) = z{27r) and initial guess 

^2(^,0) ^^i(^,ri(^)). 



Having computed ^2(^, 0), ^{9, 0), 



'',0), we obtain the approximation 



which is accepted in the region 



{{9,r) -.{] <9 < 27r, ri{9) <r <ri{9)+ es}. 



(22) 



In this way, we extend our approximation of ip{9^ r) to the annulus (l22l) and our running 
approximation is 



M0,r), 0<r<ro{9), 

ro{9) <r <n{9), 
^Md,r), n{9) < r < n{9) + €2 



for ^ e [0, 27i] 



5. Steps 4-5 can now be iterated. Indeed, suppose we have computed an approximation 
'ipk{9, r) to ip{9, r) of the form ^lJk{9, r) = \E'^ {9, r — rk-i{9)), and defined in the region 



{i9,r):0<9< 271, n-ii9) < r < rk-ii9) + ej, 
where ^k{9, a) = ^{9, rk-i{9) + a). 



(23) 



N 



^^(^,a) = vI/,(^,0) + 5^^(^,0) 



i=l 



a 



and rk_i{9) is the solution to fllSaj) with initial condition rfc_i(0) = S.^gCj. To extend 
our current approximation of i/j beyond the domain of -0^, we begin by computing the 
solution rk{9) to fllSap with initial condition rfc(O) = rfc_i(0) + e^. This is done by 
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solving a BVP using (I15ap with boundary conditions r(0) = r(27r) = rk-i{0) + and 
initial guess rfc_i(-) + e^. Let now \E'fc+i(^, a) = i/j^O, r^^O) + a) and form 

N 



'k+l\ 



i=l 



The coefficients ^-§^T^{d,0) are computed by solving BVPs using the ODEs (|2UI) with 



boundary conditions ^-^ti(0,0) 

Similarly, to compute ^^^+1(^5 0) we solve a BVP using (115bp with boundary conditions 
z{0) = z{2it) and initial guess 

Having computed \E'fc+i(6', 0), ^^^(0, 0), . . . , ^-^^{6, 0), we obtain the approximation 

which is accepted in the region 

{ie,r):0<e< 27r, rfc(^) < r < niO) + ek+i}. (24) 

To avoid overlapping the domains of ipk and "0^+1 5 the domain of definition (l23l) of ipk 
is redefined to be 

{{e,r):0<e< 27r, n-iiO) <r< rk{9)}. 

6. After iterating Steps 4-5 a A; > 1 number of times, we obtain the following piecewise 
smooth approximation to ip{0, r): 

>o(^,r), 0<r<ro{e), 



i(27r,0) and initial guesses 



(25) 



^'i/jk{9,r), rfc-i < r < rfc(6') 

for 9 E [0,27r]. 

Let us make a few remarks about our algorithm. 

Remark 3.2. 1. The coefficients |^(^,0), 2 = 0,1,..., are computed by solving BVP 
problems for mainly two reasons: (1) they are known a priori to be periodic, and (2) 
we have good approximations of them from the previously computed coefficients. 
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2. The main computational effort of our method is in computing the coefficients ^^{0, 0). 
From fl2Ul) . we see that the ODE for ^■(6',0) is linear in ^-(^,0) and is polynomial 
in the previously computed "^{9,0), ... , g'^^-? {0, 0). Hence, a way to speed up the 
computation is to solve for the coefficients ^^(^,0) order- by-order. This can result in 
computational savings when the RHS of fl20|) is complicated to evaluate or when nx N 
is large. 

3. When the tw-dynamics are given by the harmonic oscillator 

Wi = -W2 
W2 = Wi 

the computation of the radial curves r{9) is trivial and are given by the constant 
curves r{9) = r(0). In this case, there is no need to redefine the outer boundary of the 
successive approximations ipi when going from one annulus to the other. 

4. In order for the algorithm to produce a meaningful approximation to i/j, the domain on 
which the approximation (125|) is defined, namely {{9,r) : < 9 < 2it,0 <r< rk{9)}, 
must of course be contained in the cylinder 9 G [0,27r], |r| < e on which is defined. 
Since e is not known a priori, the algorithm must proceed from r^-i to r^ by taking 
small increments rfc(O) — rfc_i(0) = and choosing ro(0) sufficiently small. 

To end this section, we prove that the sequence of approximations {4'k}'^^i obtained from 
fl25|) convergence uniformly to ifj. 

Theorem 3.1. Suppose that is defined on the cylinder Q = {{9, r) : < 9 < 27r, < r < e} 
and let e < e he chosen so that if r : [0,27r] — )■ M zs a trajectory of fllSal) with r(0) < e then 
r{9) < e. Let i/jk be defined as in (125|1 . with step-size rj{0) — rj_i(0) = ej := -^e, for 
j = 0,1, . . . , k. Then ipk ^ "4^ uniformly in Q = {{9, r) : < 9 < 27v, < r < rk{9)}. 

Proof. We first note that the existence of e follows by Lyapunov stability of the exosystem. 
By construction, rj{0) = for j = 0,1,..., k, and in particular rk{0) = e, thereby 

rendering the domain Vt independent of k. Also we note that, by shrinking e if necessary, by 
Gronwall's lemma it follows that 

\r{9) - f{9) I < \r{0) - f(0) |e^^ (26) 

for all trajectories 9 r[9) and 9 f{9) of fllSap such that < r(0), f(0) < e, where K is 
a Lipschitz constant independent of 9. 

Now by definition, we have that '^j{9, a) = '>p{9, rj-i{9) + a) and therefore 
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Hence, by Taylor's theorem, 



1 f)N+^Sii . 



iV! 7o -^r^+i 

Let C = max(e ,,)ef^ II apv+T(6',T)||, and we note that C exists by continuity of q^n+i on the 
compact set fl. We therefore have that 



|^,(^,a)-vl>;(^,a)|| <C 



(iV + 1)! 



provided < cr < rj{6) — rj^i{6) for 6 G [0, 27r], for all j = 0,1, . . . , k. Now since rj{0) < e 
for j = 0,l,...,k,it follows by ([26]) that 



|r,(^) - r,_,m < |r,(0) - r,-i(0)|e^^^ < ^eV'^^. 

Therefore, given (^,r) G fi, say that rj_i(^) < r < rj(6') for some j G {0, 1, ... , A;}, it 
follows that 

r) - Md, r)\\ = \\^,i9, r - r,_i(^)) - r - r,^i{9))\\ 
= 11^,(0, r-r,_i(^))-^f(^,r-r,_i 



-(iV + l)!^^+i ^ ■ 
Hence, ||?/'(0,r) — 'i/jk{0,r)\\ — )■ as A; — > oo uniformly in Q. This completes the proof. ■ 

4 Examples 

In this section we present examples illustrating our method. 

Example 4.1. In this example we take a linear dynamical system of the form (fT2|) . whose 
center manifold is easily computed, perform a nonlinear change of coordinates and arrive 
at a nonlinear system on which we apply our method. The true solution for the nonlinear 
system is then readily available and we can compare the approximations produced by our 
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method with the true soluiton. Consider then the hnear dynamical system 



Xi = X2 + hwi + ^W2 



X2 = Xs + Iwi + ^W2 



X3 
Wi 

W2 



—Xi - 

-W2 

Wi. 



-,W2 



The center manifold equation for this system in the unknown mapping x = <^(w) is 



dw 



[w)Sw = C(/){w) + Dw 



where x = {xi,X2,X2,),w = {wi,W2), and 



S 



-1 

1 



c 




(\ 



D 



It is straightforward to verify that (j){w) = {—^Wi, —^Wi — |w2, —^wi — ^102) is the unique 
solution to the center manifold equation for this system. Consider the coordinate change 
z = Z{x) = (— 3xi, 9xi — 6x2, —X2 + X3 + p{—3xi, 9xi — 6x2)), where p : — )■ M is a smooth 
function. The system in {z, w) coordinates takes the form of (fT2|) with matrix B having 
eigenvalues —1, | ± y/Si. By direct substitution, the solution to the center manifold equation 
of the system in {z,w) coordinates is z = Z{(j){w)) = {wi,W2,p{wi,W2)), i.e., it is the 
graph of the function p. For purposes of illustration we take the egg carton shaped function 
p{wi,'W2) = sin(iyi) sin(t/;2), whose graph is shown in Figure [TJ The patchy approximation 
4'k{wi,W2) = {wi,W2,p{wi,W2)) computed with our method with k = 10 annular regions of 
thickness e = 0.5 and of order = 2 is shown in Figure [21 The error between the patchy 
approximation and the true solution is shown in Figure O Lastly, Figure S] shows that in 
order to get a similar error bound as with the patchy approximation, one needs to use a 
polynomial approximation of degree 19. 

Example 4.2. This example illustrates the loss of stability when using polynomial ap- 
proximations. We proceed as in Example 14.11 but instead use the volcano type function 
p{x, y) = sin(x^ + 2/^)6^""^^"^^, whose graph is shown in FigureO The patchy approximation 
is shown in Figure |6] and the error in using the patchy approximation is shown in Figure [3 
The patchy approximation is constructed with k = 60 annular regions of thickness e = 0.05 
and of order = 1, i.e., we only use a first order Taylor series in the radial direction. For 
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Figure 1: True solution p{wi,W2) = sin(wi) sm{w2). 





Figure 4: Error p{wi,W2) — p{wi,W2) with polynomial approximation of degree 19. 

this example, polynomial approximations of orders up to 30 where tested and it was verified 
that as one increases the order of the polynomial approximation the error in fact increases 
on the domain in consideration. 




Example 4.3. Consider the inverted pendulum cart system from Example 12.11 with two- 
dimensional exosystem given by 




where a > 0. System (p7|) is a special case of the unforced Buffing's oscillator with no 
damping [4]. The equilibrium w = of ( l27|) is a center and representative periodic solutions 
encircling w = are shown in Figure [3 
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Figure 6: Patchy solution p{wi, W2) with = 1 and k = 60. 




Figure 7: Error p{wi, W2) — p{wi, W2) with patchy solution with = 1 and k = 60. 




w 



Figure 8: Unforced Buffing's oscillator with no damping and a = j. 
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Let p{w) = —Wi. As in Example I2.H let (z,^) denote the standard normal coordinates, 
where z = (a;3,X4 + ^cos(2;3)) and ^ = {h{x), Lfh{x)) = {xi,X2). Following the notation 
at the beginning of §2], let (pi{w) = —p{w) = Wi and (p2{w) = —Lspiw) = W2- Then the 
differential equation ([8]) becomes 

Zl = Z2 - \w2 cos(2;i) 

ia = f sin(2;i) - )z2W2 sin(^i) - j^wl sin(zi) cos(2;i) 

(28) 

Wi = W2 

W2 = —wi — awf 

A patchy approximation to the solution (j){wi,W2) = {(l)i{wi,W2),4>2{wi,W2)) of the center 
manifold equation of f l28|) was computed and is illustrated in Figures iQlfTOl The data g = 10, 
^ = |, and a = J was used. The solution was computed with k = AO and radial step at the 
initial angle 9 = was taken to be = 0.05, i = 1,2, .. . ,25, and the order of the radial 
Taylor polynomials were chosen as = 2. 




Figure 9: Patchy approximation to (f)i{wi,W2)- 

The computed patchy approximation to the center manifold PDE for f l25]) is used in an 
output tracking controller of the form 

a{z, w) = k{w) + K{{z, ^) — 7r{w)) 

where 7r{w) = {(/){w),(p{w)) and k{w) = Ue{'T!-{w),w), where the gain matrix K is chosen as 
the solution to an LQR problem for the linearization of the inverted pendulum. In the LQR 
problem, the matrices Q = diag(4, 4, 4, 4) and R = 1 were chosen. A simulation is performed 
in which the pendulum is initialized at an angle of 15 degrees from the vertical and the cart 
is initialized at —0.25 from the origin. The reference trajectory, yre{{t) = wi{t), is chosen 
with initial condition j/ref(0) = 1.2. The results of the simulation are shown in Figures [mfT2l 
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Figure 10: Patchy approximation to 4i2{'Wi,W2)- 




Figure 11: Output y{t) = Xi{t) and reference Urefit) = Wi{t). 




10 15 20 25 30 35 40 
time [sec] 



Figure 12: Tracking error e(t) = y{t) — Urciit)- 



20 



5 Conclusion 



We have presented a method to compute solutions to the FBI equations of real analytic 
control-affine systems with two-dimensional exosystems. Our technique is based on the 
patchy method in [10] and on the results in [1] for uniqueness of solutions of two-dimensional 
real analytic center manifolds. In comparison with direct Taylor polynomial approximations 
pin], our method lessens the computational effort needed to produce approximate solutions 
by taking into account the periodic nature of a two-dimensional exosystem. We proved 
that our method generates a sequence of approximations converging uniformly to the true 
solution. 
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